Nonlocal imaging by conditional averaging of random reference measurements 
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We report the nonlocal imaging of an object by conditional averaging of the random exposure 
frames of a reference detector, which only sees the freely propagating field from a thermal light 
source. A bucket detector, synchronized with the reference detector, records the intensity fluctua- 
tions of an identical beam passing through the object mask. These fluctuations are sorted according 
to their values relative to the mean, then the reference data in the corresponding time-bins for 
a given fluctuation range are averaged, to produce either positive or negative images. Since no 
correlation calculations are involved, this correspondence imaging technique challenges our former 
interpretations of "ghost" imaging. Compared with conventional correlation imaging or compressed 
sensing schemes, both the number of exposures and computation time are greatly reduced, while 
the visibility is much improved. A simple statistical model is presented to explain the phenomenon. 

PACS numbers: 42.50.Ar, 42.30.Va, 42.50.St 



Classical image formation is most commonly realized 
by recording the interaction information between a radi- 
ation source and the object onto a detector, all along a 
single light path. However, in the technique now known 
as "ghost" imaging (GI), two beams are used from the 
same optical source; one beam interacts with the target 
and is collected by a so-called "bucket" detector which 
has no spatial resolution and only measures the total in- 
tensity, while the other propagates through free space to 
a reference detector which does have spatial resolution. 
Neither detector can "see" the object on its own, but 
when their second-order correlation is measured, the im- 
age appears in the coincidence data. This seems quite 
contrary to intuition, and after the first experiment was 
demonstrated with entangled twin beams [1], the term 
"ghost" was coined to emphasize its peculiar nonlocal 
nature, which was attributed to the quantum properties 
of the biphoton source. Because of the potential applica- 
tions, interest soon spread to other sources [2[, followed 
by a profusion of studies on GI with pseudothermal light 
p-jl, true thermal light [H, and even GI without a lens 
PJJ-[l5j. The expression "nonlocal imaging" is now un- 
derstood to mean that the light recorded by the refer- 
ence detector travels through free space and never inter- 
acts with the object, and is still generally accepted in 
the GI community irrespective of whether the source is 
quantum or classical. Since the prerequisite for corre- 
lated imaging is the identical spatial distribution of the 
field intensities in the two beams, it was then realized 
that the reference beam could just as well be replaced 
by an artificially generated random spatial modulation 
of the object beam, thus only a bucket detector would 
be necessary. After the experimental demonstration of 
this kind of computational ghost imaging [l6|, [I?) , atten- 
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FIG. 1: Schematic of correlated imaging. BS: beamsplitter, 
(a) Conventional GI; (b) Correspondence imaging with se- 
lected reference information. 



tion turned to improving the techniques involved in the 
data collection, correlation and averaging of thousands 
of exposure frames. Compressed sensing was found to be 
very effective in greatly reducing the number of frames 
required when the image data is sparse, and has already 
been implemented [l8[ . 

Despite all these advances, there is still much debate 
on the various quantum and classical interpretations of 
GI 0, G3 • Basically, two different interpretations have 
been offered: one is multi-photon interference based on 
Glauber's quantum optical coherence, and the other is in- 
tensity fluctuation correlation. From the point of view of 
image processing, the key to GI reconstruction is the ex- 
traction of information according to the time correspon- 
dence between the bucket signals and reference frames. 
From the mathematical viewpoint, the image retrieval is 
akin to an inverse problem in which the features (pix- 
els) of an object are to be inferred, given a time series of 
bucket and reference light field intensities. 

The basic scheme of a correlated GI system is shown in 
Fig. [Ha). Thermal (or entangled) light is separated into 
two parts through a beamsplitter (BS). An object char- 
acterized by the transmission function T(a?o), where xo 
indicates the transverse coordinate in the plane of the ob- 
ject, is inserted in front of a bucket detector Db (without 
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any spatial resolution), which registers the total intensity 
1b transmitted (or reflected) by the object, while a ref- 
erence detector Dr with spatial resolution is in the other 
empty arm. The distances from the source to the object 
and to Dr are ds and dR, respectively. For lensless GI 
with thermal light, the phase-conjugate mirror require- 
ment cIb = cIr must be satisfied [71 llQj| . In conventional 
GI, an image of the mask is retrieved from the coincidence 
measurements of Dr and Dr when the latter is scanned 
in the corresponding image plane, so all the information 
from both the detectors is required. The basis of GI is 
the transverse spatial second-order correlation [20| 
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g^{xbA) = 1 Yl 

{ti\i=l,--- ,N} 



Ib (xb,U)Ir (xR,t$l) 



Here Ib (xb,U), Ir (xr , U ) and xr,xr are the intensities 
at time ti and the transverse coordinates at Dr and Dr, 
respectively; N is the number of the coincidence measure- 
ments, i.e. exposure frames. The output Ib(U) of the 
bucket detector is the spatial integral J dxslB (xb,U), 
which erases the spatial resolution of the target. 

Below, we demonstrate a method using incoherent 
thermal light by which an image of the object can be ob- 
tained very simply from the raw reference detector data, 
without the need of any complicated image reconstruc- 
tion method such as that based on compressed sensing 
[21] . In this method, which we shall call correspondence 
imaging (CI), we only have to perform signal averaging 
over those reference frames selected by an appropriate 
gate from the bucket detector. Since no direct second- 
order correlation is involved, the theory seems to chal- 
lenge all our previous interpretations of GI. In addition 
to the normal positive image, a reversed negative image 
can also be obtained, and the quality of the images may 
be far superior to that of conventional GI for the same 
number of frames. 

A schematic of the CI setup is shown in Fig. QJb). As 
in (a), a linearly polarized He-Ne laser beam diameter 
5 mm is projected onto a rotating ground-glass disk to 
produce pseudothermal light [22|, which is then split by 
a 50:50 non-polarizing BS into two beams, one of which 
passes through an object mask to detector Db, and the 
other goes directly to Dr, to produce the bucket inten- 
sity Ir and the reference intensity distribution Ir (xr), 
respectively. The distance from the ground-glass disk to 
BS is 19 cm, and from BS to the mask 11 cm, that is, 
dp = d,R ~ 30 cm. Both detectors (Imaging Source DMK 
31BU03) are synchronized by the same trigger pulse at a 
rate of about 3 Hz. The data was acquired with an expo- 
sure time (0.1 ms) much shorter than the coherence time 
of the laser, and saved to a computer. A total of 50, 000 
frames was grabbed by each camera per image plot. The 
bucket detector could also be just a point detector so long 
as it collects all the light transmitted through the mask, 
for example by means of a focusing lens. 

The integrated beam intensity Ib(U) a t the bucket de- 
tector recorded as a function of time is shown in Fig. [2(a); 
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FIG. 2: Bucket detector intensity Ib- (a) Integrated beam 
intensity Ib recorded as a function of time sequence, (b) 
Probability distribution of Is centered about its mean. The 
left (blue) block is binned for AIb < 0, and the right (red) 
block for AIb > 0. The dashed line denotes Ib- (c) Multi- 
section binning. From left to right, the dashed lines denote 
In, Ib and ip, where In and Ip satisfy P(In) = P(Ip) — 
P{Ib)/2. 



it can be seen that the intensity fluctuations are ran- 
dom. We then calculate the deviation from the mean 
AIb(U) — Ib(U) — Ib, and plot the histogram of Ib(U)i 
or equivalently of the normalized intensity fluctuations 
AIb(U)/Ib, to show the probability distribution P(Ip) 
of Ib- Using the intensity Ib(U) a t time U as an indica- 
tor, we divide {ti} into two subsets: 

t- = {U \I B (U) <I B },t+ = {U \I B (U) >I B }- (2) 

All the Ip(U) with U G t- (or t + ) contribute to the left 
(or right) half of Fig. El(b). 

Because both detectors are synchronized in time, the 
reference signals Ir (xr, ti) may also be divided into two 
groups according to U G or U G L. This is equiva- 
lent to labeling Ir (xr, U) with the time stamp of Ip(U). 
By dividing the range of Ib into more sections, say at 
In and Ip as well as at Ip, with AIn < 0, Alp > 0, 
and P(In) = P{Ip) — -P(7b)/2, we can perform multi- 
section binning for Ip (Fig. [2(c)), and hence, correspond- 
ingly, for I R (x r ,U). 

Due to the free propagation of the thermal light field, 
each reference frame is purely random, so the average 



value of all the reference signals is 

Ir (xr,U) 



Ir (x r ) = 1 

{U\i=i,- 



Ir. 



(3) 



Statistically, the mean value of a random variable is 
a constant, so we would not expect Ir (xr) to reveal 
anything, however long the exposure time, as shown in 
Fig.[3ja), and for a homogeneous thermal field each pixel 
would have the same averaged intensity. But amazingly, 
after the above partitioning, the target in the object arm 
can be reconstructed merely by taking the average of only 
the reference intensities within a selected distribution: 

1 ^ Ir(x r ,U), (4) 



R±(x R ) = — £ 

{ti\uet±} 



where N± = ^2 t . et± 1 is the cardinality of t± (see 
Figs.^b) and (c))* 

We know from information theory and statistics that 
a complete random basis is recommended for image re- 
construction. In our experiment this is supplied by the 
spatial fluctuations of all the random modes of the ther- 
mal beams. The object and Dr are at the same distance 
from the source, so they see the same mode distribution. 
Part of this field, i.e. the transmitted (or reflected) light 
through the mask, contributes to the bucket signal, and 
may fluctuate above or below its mean value in each ex- 
posure. The intensity of the corresponding Dr pixels 
will also fluctuate in harmony, but the rest of the frame 
will be fluctuating in an uncorrelated way. Let us divide 
the spatial pixels at the mask {xo} into subsets X\ and 
Xq with the former transmitting and the latter blocking 
light. It is the pixels {xo \ xo E Xi} that contribute 
to the bucket signal Ib , and when the field intensity 
fluctuates higher (lower) in this area, the bucket output 
will naturally increase (decrease). The intensity of the 
corresponding reference detector pixels will also increase 
(decrease), but as this contribution is superimposed upon 
the entire beam intensity, i.e. all the light from both ar- 
eas Xi and Xo, the total reference output will be more 
or less constant. By separately collecting and averaging 
over a sufficiently large number of exposures according to 
the positive (negative) bucket signals, the positive (neg- 
ative) image of the object {X\) will then stand out from 
the constant background. 

A conceptual mathematical model is as follows. Since 
the field distribution at the mask is the same as that at 
Dr, the bucket intensity may be written 



Ib (U) = J2{Z}Ir (x, U) \T (a 



(5) 



where for convenience we have dropped the subscripts 
of x. From Eq. © we thus have I B = IrY, I t (^)| 2 - 

{x} 

For a narrow bucket intensity bin, centered at some 
Ib{U) 

approximation 



lg far above we may make the following 



Ib(U) ~ Yl Ir&U)\T(x) 

{x\xeX!} 



(6) 
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FIG. 3: Images obtained by nonlocal CI; all figures are au- 
tomatically gray-scale compensated, (a) Average of all the 
frames from Dr. (b) Negative and (c) positive correspon- 
dence images after time binning of Dr based on negative and 
positive AIb, respectively, (d) Images obtained using all the 
frames from Dr but with the negative image signals R- (x) 
inverted, (e) and (f), as in (b) and (c) above, respectively, 
but after normalized array division by (a). 



Since Nb+ is the number of Ib{U) values registered in 
the range of we can write 



G+(x) 



N B + 



Ib(U)I r (S,U) 

{U\I B (U)^I+} 



N B + 



{U\I B (U)^I+} 



A 2 



On the other hand, by using Eq. (j5j), we have 

G+(x) ^ I b Ir 

where A 2 R = £ 

{U\i b (u)^i+} 
means that, for x G X\ 



r\T(x)\\ (8) 
[I R (x,U)-I R ] 2 . This 



I+(x) ocC b + \T(i 



(9) 



with Cb = IbIr/^r being a constant due to the back- 
ground. Similarly, for x G Xo we have I+(x) oc Cb + 0. 
For a narrow bin centered at some 1% far below 
the same argument is applicable under the exchange 
X\ <H> Xo. Thus, while the summation of reference 
frames corresponding to Ib(U) ~ Ib gi yes the posi- 
tive image of pixels in Xi, the summation frames for 
Ib(U) ~ Ib gi yes the negative image of pixels in Xo- 

When going from the narrow bin at to a wide bin 
of Ib (U) > Ib 5 errors from two sources come in: since 
Ib (U) ^ Ib is approximated by the single value i^, the 
validity of Eq. (j6j) becomes weaker. This implies that the 
reference frames with their corresponding Ib (U) close to 
Ib contain limited information about |T (x)| 2 , relative to 
the background, so the image quality is reduced. How- 
ever, if we take R + {x) — R-(x) then the background can 
be removed. We may also define the normalized distri- 
butions Jzp(x) = R t (x)/Ir(x). 

Figure [3] shows the images obtained for the letters 
"CAS" in an object mask which is 1.4 x 2.56 mm 2 in size 
(300 x 550 pixels). We define two levels of transmission, 
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1 for total transmission through the letters, and else- 
where. The results obtained by various different methods 
are presented for comparison to illustrate the advantages 
of CI. We see from Fig. [3]^a) that, as expected, merely 
taking the time average Ir (x) using all the 50, 000 frames 
of the reference detector produces nothing (the slanted 
lines are artifacts from the rotating ground glass plate 
due to the long time exposures). After partitioning by 
the two complementary time series given by Eq. (j2j) , there 
emerges a negative image i?_ (x) from 26, 005 reference 
frames (Fig. |3fb)), and a positive image R+ (x) from the 
remaining frames (Fig.[3jc)). The image R+ (x) — R- (x) 
obtained using all the information from Dr but with the 
negative image signals inverted is shown in Fig. [3fd). 
When we take the normalized correspondence images 
from Dr, obtained as 7_ (x) and 7+(x) by using ma- 
trix array division, we obtain the images in (e) and (f), 
respectively, which are much clearer than in (b) and (c) 
above. It is interesting to note that in all the upper row 
figures, there is a small black dot with a white centre 
at the bottom of the letter "C", due to diffraction from 
a dust particle on the surface of the reference detector. 
However, this is no longer visible in the lower row as the 
background has been removed through normalization. 

It should be noted that although CI still requires syn- 
chronization with the bucket detector signals, it is quite 
different from the usual correlated imaging of GI. Each 
individual frame of Dr is random and the object cannot 
be seen from the mean of the total, but its image can be 
retrieved from the conditional average of subsets of the 
data, a seeming contradiction. Moreover, the visibility 
of the images may be optimized by appropriate selec- 
tion and weighting of the data partitioning, for example, 
through appropriate choice of bucket detector intensities 
In and ip, we can divide the histogram into four areas, 
as illustrated in Fig. [2^c). Besides 7p, we introduce two 
more In (< Ib) and Ip (> Ib) to divide the range of 
Ib into four subintervals: Pi = [0,/at), Pi = [In,Ib), 
P3 = [Ib,Ip)i and P4 = [Ip, +00). Assuming that there 
are Ni measurements (or frames) with a bucket intensity 
within the range Pi (i = 1, 2, 3, 4), in the same way as we 
defined R T (x), we introduce 

Ri{x) = h&U)' ( 10 ) 

* {u\i B (U)ePi} 

and its normalized form ji(x) = R%{x)/ Ir(x). We expect 
that a negative image would be obtained with R\ or R2 
and a positive image with R3 or P4. 

Using the same experimental data as above, three sets 
of values for Ni were chosen to create the correspondence 
images of Fig. HI In the first row, the number of refer- 
ence frames contained in each sector Pi was the same 
(Ni ~ N2 ~ N3 ~ 7V4). It is obvious that the images 
obtained from Pi and P4 are far superior in quality com- 
pared to those from P2 and P3, although the same num- 
ber of frames were used. In the second row, a weighting 
of TVi : 7V 2 : ^3 : ^4 ~ 1 : 9 : 9 : 1 was used. Even 
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FIG. 4: Selected CI images when the information from Dr is 
divided into four subsets with different weightings (the gray- 
scale has been automatically adjusted). (I) Each subset is 
composed of the same number of measurements, Ni ~ iV~2 ~ 
N3 rsj N4; (II) most of the frames are concentrated in the 
central areas P2 and P3, with a weighting of Ni : N2 : N3 : 
iV 4 ~ 1 : 9 : 9 : 1; (III) when JVi : N 2 : N 3 : N 4 ~ 4 : 1 : 1 : 4, 
no images can be deciphered from P2 and P3. 



in this case, we see that the reconstructions from the ex- 
terna, i.e. the wings of the fluctuation distributions in 
Figs. Efb) and (c), are comparable with those from the 
central area, despite the fact that only 5% of the data 
was used. In the third row, the ratios of the four parts 
are about Ni : N 2 : N 3 : iV 4 ~ 4 : 1 : 1 : 4, and we 
can no longer see any image in the regions of P2 and 
P3. This indicates that most of the object information is 
concentrated at the two ends P\ and P4, i.e., the refer- 
ence frames corresponding to larger bucket intensity fluc- 
tuations contain more information and contribute more 
to image retrieval. When the image is indistinguishable 
from the background, we may infer that the exposure 
intensities resulted from both subsets X\ and X$. By re- 
fining our weighting of the partition selection, the image 
may be recovered faster and with fewer measurements, in 
agreement with the experimental observations of Fig. 2J 

It should be noted that in Figs. [3] and |4] the gray-scale 
has been automatically adjusted by the software for bet- 
ter pictorial detail: the minimum of the image matrix is 
displayed as the darkest and the maximum as the bright- 
est so that contrast is artificially improved. In this way, 
images that originally had very different visibilities may 
look quite similar on a computer screen, or vice versa. 
The CI images appear comparable to those obtained by 
conventional GI although only half the data from Dr or 
even less is used. 

To demonstrate the superior visibility of CI more 
clearly we also performed an experiment under the same 
conditions as before but with a simple double-slit as ob- 
ject (slit width and distance 0.2 mm and 0.4 mm, re- 
spectively). The fluctuation histogram of Ib was divided 
into four sectors, and to spotlight the differences in vis- 
ibility we look at the intensity profiles along one trans- 
verse dimension. For a better comparison of the original 
image contrast of the various techniques, in Fig. [5] we 
have adopted a common gray-scale in which the image 
matrices are divided by their maximum value and dis- 
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x (pixel) x (pixel) 

FIG. 5: Image intensity profiles for a double-slit object; trans- 
verse coordinates are the pixels of Dr (each pixel is 4.65 by 
4.65 /im 2 ). Each plot is normalized by its maximum value, 
(a) Normalized negative correspondence images reconstructed 
for 7- and two 71 values, (b) Top curve: second-order nor- 
malized GI; lower three curves: normalized correspondence 
images reconstructed for 7+ and two 74 values. 

played with zero as the darkest and 1 as the brightest. 
Fig. [5fa) shows the normalized negative correspondence 
images, where 7_ indicates the negative image when the 
probability histogram is only divided into two parts, and 
the percentage values in 71 (25%) and 71 (5%) represent 
the percentage N\ of reference frames selected in parti- 
tion Pi. Fig. Efb) shows positive images: the top curve 
is from normalized second-order correlation GI; 7+ is the 
normalized positive CI image, and 74(25%) and 74(5%) 
are the CI plots for N A /N ~ 25% and 5%. The superi- 
ority of the latter is immediately obvious. Moreover, the 
less the number of exposures that are taken, the better 
is the visibility relative to that in conventional GI. The 
definition and analysis of the visibility and signal-to-noise 
ratio in CI is a complicated problem [23|, [24| , so a more 



quantitative discussion, as well as a full theory of CI in 
general, will be presented in a future paper. 

In summary, we have demonstrated the reconstruc- 
tion of positive and negative images of a nonlocal object 
through selective averaging of the exposures of a reference 
detector that has never interacted with the target field. 
A simple statistical model is proposed to explain the phe- 
nomenon, which, to our knowledge, cannot be explained 
by conventional classical or quantum wave optics. Since 
the method is based on a form of conditional averaging 
of the reference data, much less information is required 
than in GI, which in fact contains much redundant infor- 
mation. Moreover, the visibility is significantly better, 
especially if the number of exposure frames is limited. 
The reconstruction process only involves straightforward 
selection and stacking, so it is much simpler than conven- 
tional GI or compressed sensing methods; complicated 
calculations and algorithms are not required, and stor- 
age space, memory consumption and processing time are 
greatly reduced, which is a particular advantage when 
the images are large. The basic concept may also be 
extended to similar experiments that use correlation cal- 
culation, including computational GI, so potential ap- 
plications are plentiful. Although this CI technique is 
apparently simple, the underlying physics is quite subtle 
and deserves further exploration. It could open up new 
opportunities in the field of metrology, positioning, and 
imaging. 

This work was supported by the National Natural Sci- 
ence Foundation of China Grant No. 60978002, the 
National Program for Basic Research in China Grant 
Nos. 2010CB922904 and 2007CB814800, and the Na- 
tional High Technology R&D Program of China Grant 
No. 2011AA120102. 



[1] Pittman T B, Shih Y H, Strekalov D V and Sergienko A 

V 1995 Phys. Rev. A 52 3429(R) 
[2] Cheng J and Han S S 2004 Phys. Rev. Lett 92 093903 
[3] Bennink R S, Bentley S J and Boyd R W 2002 Phys. Rev. 

Lett. 89, 011360; Bennink R S, Bentley S J, Boyd R W 

and Howell J C 2004 Phys. Rev. Lett. 92 033601 
[4] Gatti A, Brambilla E, Bache M and Lugiato L A 2004 

Phys. Rev. Lett. 93 093602 
[5] Ferri F, Magatti D, Gatti A, Bache M, Brambilla E and 

Lugiato L A 2005 Phys. Rev. Lett. 94 183602 
[6] Cai Y J and Zhu S Y 2005 Phys. Rev. E 71 056607 
[7] Cao D Z, Xiong J and Wang K G 2005 Phys. Rev. A 71 

013801 

[8] Valencia A, Scarcelli G, DAngelo M and Shih Y 2005 

Phys. Rev. Lett. 94 063601 
[9] Zhang D, Zhai Y H, Wu L A and Chen X H 2005 Opt. 

Lett. 30 2354 

[10] Scarcelli G, Berardi V and Shih Y H 2006 Appl. Phys. 

Lett. 88 061106 
[11] Basano L and Ottonello P 2006 Appl. Phys. Lett. 89 



091109 

[12] Meyers R, Deacon K S and Shih Y H 2008 Phys. Rev. A 

77 041801(R) 
[13] Liu H L and Han S S 2008 Opt Lett. 33 824 
[14] Ferri F, Magatti D, Sala V G and Gatti A 2008 Appl. 

Phys. Lett. 92 261109 
[15] Chen X H, Liu Q, Luo K H and Wu L A 2009 Opt. Lett. 

34 695 

[16] Shapiro J H 2008 Phys. Rev. A 78 061802(R); Erkmen B 

1 and Shapiro J H 2010 Advances in Optics and Photonics 

2 405 

[17] Bromberg Y, Katz O and Silberberg Y 2009 Phys. Rev. 
A 79 053840 

[18] Katz O, Bromberg Y and Silberberg Y 2009 Appl. Phys. 

Lett. 95 131110 
[19] Shih Y 2007 IEEE J. Sel. Top. Quantum Electronics 13 

1016 

[20] Glauber R J 1963 Phys. Rev. 131 2766 

[21] Wu L A and Luo K H 2011 AIP Conf. Proc. 1384 223 

[22] Goodman J W Statistical Optics (Wiley Classic Library, 



6 



New York, 2000); Goodman J W Speckle Phenomena in 
Optics: Theory and Applications (Roberts & Company 
(Englewood, Colorado), 2007) 
[23] Basano L and Ottonello 2007 Appl. Opt 46 6291 



[24] Brida G, Chekhova M V, Fornaro G A, Genovese M, 
Lopaeva E D and Ruo Berchera I 2011 Phys. Rev. A 83 
063807 



